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1. Introduction 

The problem of two interacting electrons in atoms, ions and molecules is of fundamental 
importance in quantum mechanics, the understanding of the Helium atom being one of 
the prime successes of the "new quantum mechanics" [l]. If the Coulombic interaction 
between electrons and nucleus is replaced by a harmonic confining potential, the system 
is referred to as Harmonium, quantum dot Helium or Hooke's law atom. After early 
work without additional magnetic field [2], this problem, with and without a magnetic 
field, has again come in the focus of interest in the past 20 years mainly for two different 
reasons: 

Firstly it is of principal interest to find (exact) analytic solutions for a problem 
with genuine Coulomb interaction. That such exact solutions do exist has been shown 
for two electrons in a harmonic confinement in 3 dimensions, without an additional 
magnetic field [21 13], and in 2 dimensions with a magnetic field by Vergin |1] and Taut 
[5]. The existence of closed solutions for specific ratios of the Coulomb interaction 
and the harmonic potential is related to an accidental (hidden) symmetry visible after 
mapping the Coulomb interaction to a four dimensional harmonic oscillator potential 
[6]. Furthermore, in 2 dimensions and using perturbation theory for the Coulomb 
interaction, analytic solutions have been given in [7]. The analytic solutions can be 
put to good use in the testing of density functionals and other methods for many-body 
systems [H [9l [10]. Secondly, due to recent progress in nanofabrication, few electron 
quantum dots have come into the limelight as they may provide a realization of a 
quantum bit [TT] . 

Due to the fact that exact analytic solutions of the Schodinger equation are 
limited to an infinite set of discrete oscillator frequencies [Sj [5], also approximate 
analytic solutions have been sought for. Apart from the perturbation theoretic ones 
in 2d, mentioned above, also semiclassical approaches have been taken which rely on 
the Wentzel-Kramers-Brillouin (WKB) |12], respectively the Einstein-Brillouin-Keller 
(EBK) quantization rules |13]. Most notably, in 2d and with an additional magnetic 
field this has been done analytically by Klama and Mishchenko [H]. It turned out 
that the problem of two electrons in 2d, although one of 4 degrees of freedom (DOF), 
is highly separable and therefore in the end only a 1 DOF WKB approach is needed, 
which describes the spectrum with a high accuracy. In another semiclassical analysis 
of the harmonium problem, including a magnetic field, dimensional effects have been 
discussed by Nazmitdinov and collaborators [T5] . 

Whereas in WKB energy information is directly available, using the semiclassical 
initial value representation (IVR) of the quantum mechanical propagator of Herman and 
Kluk (HK) [161 [13 dH] spectra can be gained by Fourier transformation of a time series, 
usually an auto-correlation function [19]. In [20] it has been shown that accurate spectra 
for coUinear Helium can be obtained in this way. Furthermore, semiclassical work on one 
electron systems in external laser fields gives the explanation for the plateau formation 
in high harmonic spectra [21] . 
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In this paper we will show that the time-dependent semiclassical approach, 
complementary to the energy-domain WKB method, is also capable of reproducing 
high quality spectra for 2d Harmonium in a homogeneous magnetic field. Thus adding 
to the knowledge of semiclassics for interacting many particle dynamics is not a mere 
exercise due to the fact that the time-dependent approach is not principally restricted 
in dimensionality nor to the case of a circular quantum dot and also additional time- 
dependent potentials may be treated in a similar fashion without much more effort. 
In Sec. O the Hamiltonian is briefly introduced, whereas in Sec. [3] we review the 
semiclassical IVR of Herman and Kluk based on frozen Gaussian wavepackets. The new 
results we gained are compared with full quantum as well as with (Langer-corrected) 
WKB results in Sec. HJand we give conclusions and an outlook in Sec. [51 



2. The Hamiltonian 



The Hamiltonian for the problem of two charged particles (charge q) in a 2d circular 
quantum dot (dielectric constant e) with a confinement frequency of Uq inside a magnetic 
field, derived from a vector potential A(rj) is given by 



2m* 

i — J. \ / I - — I 

Here k = g^/(47reeo) and for two electrons m* is the (effective) mass and the charge 
is g = — e. We are not explicitly taking into account electron spin, thereby neglecting 
the Zeeman energy of the spins in the magnetic field, see also [15]. In addition, as 
will be seen and commented on below (see Section 4.2), for the propagation in time, 
we will consider an unsymmetrized position-dependent part of the quantum mechanical 
wavefunction. 

To make progress, center of mass (cm) and relative coordinates are introduced, 
according to 

R=^(ri + r2) P = (pi + P2), (2) 

respectively 

r = ri-r2 p = ^(pi-p2). (3) 

In their terms, the Hamiltonian is separable (i. e. depends on the sum of two terms 
which depend only on cm or relative coordinates, respectively) and for A(rj) depending 
linearly on coordinate reads 

H = Hem + Hj-cl (4) 



with 



Hem ^ (P + 2eA(R))2 + ^MulR' (5) 
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where M = 2m* is the sum of the two masses and /i = m*/2 is the reduced mass. 
A constant magnetic field in the z direction is derivable from a vector potential 

A, = -yB/2, Ay = xB/2, A, = (7) 

given in Cartesian coordinates and in symmetric gauge. For the relative motion 
Hamiltonian in 2d this leads to 

^rei(p, r) = ^ {pI + pI) + 2/^^'r2 + ujlL, + -, (8) 

where Lz = pyX — ypx is the 2;-component of the angular momentum, ujl = eB/{2m*) 
is the Larmor frequency and = uJq + l>j\. 

In the following we will exclusively concentrate on the relative motion, due to 
the fact that the solution of the cm problem is trivial (it can easily be gained from the 
solution of the relative motion by neglecting the Coulomb term and a suitable rescaling). 



3. A semiclassical initial value propagator 

Following early work by Heller [16] , who first presented a semiclassical integral expression 
for a time-evolved wavefunction in terms of fixed width (frozen) Gaussian wavepackets, 
Herman and Kluk have derived the semiclassically correct prefactor for an A^-degree of 
freedom system [1^ . The resulting expression for the propagator in the case A^ = 2 is 

ir«^(r, t- r', 0) = / v/d^(r|^(p„ q,)) 

exp{l5(p',q',t)}(^?(p',q')|r'), (9) 

where 

Siv'. q', t) = f dt'[pTq - i7rei(p, q)] (10) 
Jo 

is the classical action. With a series of papers in 1994 [HI |22l [23], Kay has laid the 
ground for a fiurry of publications using the Herman-Kluk propagator. One of its main 
features is the integration over multiple Gaussian wave functions 



(r|^(p,q)) = 




with constant width parameter matrix 7 (therefore the term "frozen" is frequently used), 
which we here assume to be the same for the initial and the final Gaussian. 

Furthermore, the matrix h, whose determinant appears in the preexponential factor 
of the semiclassical propagator is given by 

h(p', q', t) = ^(mii + 7m227"^ - «%m2i - ^mi27~^) (12) 

for diagonal 7-matrix, and contains sub-blocks of the so-called stability (or monodromy) 
matrix to be discussed in more detail below. In a numerical implementation, the branch 
of the square root in Eq. (Q has to be chosen in such a fashion that the whole expression 
is a continuous function of time [T8] . 
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Classical dynamics enters the HK IVR of the propagator via the classical 
trajectories (p^ = p{p',<i',t),q^ = q(p',q',t)), that are initial value solutions of 
Hamilton's equations, and which in the case of the relative motion of 2 electrons in 
a quantum dot read 

qx = -Px - uJiqy (13) 

1 , , 

% = -Py + ^iqx (14) 

Px= - ^LPy - P^^qx + (15) 
Py = ^LPx - fJ'^^qy + (16) 

An in depth study of the classical dynamics of two interacting particles in a magnetic 
field in two dimensions has been given by Curilef and Claro [2l]. Even without the 
harmonic confining potential, and in spite of the presence of the repulsive Coulomb 
interaction, the motion is bound due to the presence of the magnetic field. 

Furthermore, also the equations of motion for the stability matrix are classical 
equations. They will be given in vector form and for = 2 we define = 
(5p^, 5p'y), ^q' = {^q'xi Sq'y) as small initial deviations in phase space. Their time evolved 
counter parts are 

= mii5p' + mi2(5q' (17) 
(Jq^ = m2i5p' + m225q', (18) 

where the 2x2 matrices mij{i,j = 1,2) are submatrices of the stability (or monodromy) 
matrix 

/ _ _ \ / ^ ^ \ 

(19) 

The equation of motion for M can be obtained by linearizing Hamilton's equations for 
the deviations and reads (for more details, see, e.g., App. 2C in [25] ) 

M = I a5 S5 I M = -JHM. (20) 



M = i 


^ mil 


mi2 ^ 






V m2i 


m22 y 






dt 



dPtdPi wqT 

Here the skew symmetric matrix 



(21) 



and the Hessian, i. e. the matrix containing the second derivatives of the Hamiltonian 

^rcl(Pi,qi), 



H-l "PS; I (22) 
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have been used. We note that the Hessian of the Hamiltonian in contrast to the 
standard case oi H = T(p) + V{q), for ^ also contains nonzero mixed second 
derivatives. 

The initial conditions follow from the definition of the stability matrix to be 



and in the numerics we solve the stability equations along with the trajectories by using 
a symplectic integrator of second order, the so-called symplectic leap frog (or position 
Verlet) algorithm [26]. This automatically ensures that the action is discretized using a 
mid-point rule, which in general is necessary in the presence of a vector potential |27] . 

The HK propagator is determined entirely with the solution of classical initial 
value problems, hence it is frequently referred to as an initial value representation, in 
contrast to the well-known van Vleck-Gutzwiller propagator [281 EH] that is based on 
the solutions of a classical root search problem (double sided boundary value problem). 
In most numerical applications published so far, the bare propagator has been used to 
evolve a Gaussian wavefunction in time (see also below). Applying the propagator to 
such a wavefunction, centered around (pQ,,qQ,) and with the same width parameter as 
the frozen Gaussians (fTTj) used for the propagator, according to 



is eased by the fact that the integral over r' can be done analytically by using 



for the overlap between the initial Gaussian wavefunction and the frozen Gaussian. 
After inserting the propagator together with Eq. fl2^ into Eq. fl2^ . an integration 
over initial phase space (the initial conditions of the classical trajectories) remains to be 
done. The sampling of the initial phase space is done using the Monte-Carlo integration 
technique [30] and facilitated by the exponential damping of the integrand far away from 
the center of the initial wavefunction. In the numerical results to be presented below 
the convergence of the time signals was checked by increasing the number of propagated 
trajectories. For four phase space degrees of freedom typically 10^ trajectories are needed 
for converged results. In cases of strongly chaotic dynamics (not considered here) and 
if long-time information is needed this number may, however, increase dramatically 
[3T] . We note in passing that in contrast to the multi-trajectory Herman-Kluk method 
[32] . the single trajectory Thawed Gaussian Wavepacket method [33] uses just a single 
trajectory to express the final wavefunction. 




(23) 




(24) 



(^7(p',q')l*«(0)) = exp{--(q' - qj^7(q' " qj 

+ 7i-(q'-qa)^(p' + pJ 




(25) 
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4. Semiclassical spectra 

In the following we compare spectra that have been gained by two different methods, 
the energy-domain WKB method and the time-domain IVR method with subsequent 
Fourier transform, to full quantum mechanical ones. 

4.1. WKB spectra 

In the case of a circular quantum dot in two dimensions that we consider, the 
semiclassical WKB approach becomes particularly attractive, due to the fact that also 
the relative motion is separable in polar coordinates |14] . 

In the simple case without the Coulomb term, a one DOF WKB quantization of 
the radial motion leads to the Fock-Darwin energies [SI [351 El US] 

Ere\{nr,m) = (271.^ + \m\ + l)hQ — mhuL (26) 

where 

= 0,1,2,... m = 0,±1,±2, ... (27) 

are the radial, respectively azimuthal quantum numbers. We note in passing that these 
eigenvalues have originally been gained for a single electron in a magnetic field and 
an additional harmonic confinement. In the case Uq = 0, i. e. without the external 
confinement, for m positive or zero, the spectrum reduces to the so-called Landau 
levels Ej-ei = {2nr + 1)Tiujl, which are infinitely degenerate with respect to the angular 
momentum quantum number, see also |12j . 

Separability is still given, also with the Coulomb term. The WKB energies can then 
be gained by solving for the roots of a quartic equation and numerically integrating a 
complete elliptic integral [Ml [15]. We stress that in the course of the derivation the 
two-dimensional Langer modification [36] of the azimuthal (magnetic) quantum number 
— 1/4 — )■ has to be performed. 

Graphical representations of the eigenvalues as a function of magnetic field strength 
can, e. g., be found in [71 [TH [15]. It is worthwhile to note that as the magnetic field 
increases, the ground state shifts to levels with higher angular momentum because the 
(repulsive) Coulomb energy gets smaller when the angular momentum grows along with 
the average distance between the electrons [H]. Here, we give a collection of some 
eigenvalues in table [H where by looking at the results with the Coulomb potential, one 
can see that for the magnetic field chosen, the m = 1 state is now the ground state. For 
the numerics in this paper, for reasons of simplicity, we made all quantities dimensionless 
by setting /i = ^ = e = /€ = l (please note that this does not correspond to atomic 
units). 
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without Coulomb 


with Coulomb 




m = 


m = 1 


m = 2 


m = 


m = 1 


m = 2 





1.41 


1.83 


2.24 


3.14 


2.84 


3.02 


1 


4.24 


4.66 


5.07 


5.75 


5.55 


5.78 



Table 1. (Dimensionless) Fock-Darwin WKB energies (k = 0) and WKB energies 
with the Coulomb repulsion (k = 1) for lul = loq = 1. 



4-2. Spectra from time-series 

Now we turn to the determination of energies from time series. By using the semiclassical 
propagator an auto-correlation function can be calculated, according to 

c{t) = {mmt))- (28) 

This time series can then be Fourier-transformed into the energy domain. Using 
the expansion of the wavefunctions in (orthogonalized) energy eigenstates |^I/(t)) = 
'EnCn\n) exp{-iEnt/h}, we get 

S{cu) = / dte'"*c(t) (29) 



27Th 

oo 

Y,\cr,]^5{En-hcu). (30) 



n=0 



The peaks of the spectrum S{u) are thus located at the eigenvalues of the Hamiltonian 
[IHl [25]. In principal also a nonlinear procedure called harmonic inversion can be used 
to extract the spectral information. This has the advantage that shorter time signals 
and therefore less trajectories in the semiclassical calculations are needed [37] . 

We have performed the procedure just outlined both semiclassically with the HK 
propagator as well as fully quantum mechanically, using a finite difference method [38] 
and Cartesian coordinates for the two DOF. A comparison of the two different time 
series for an initial Gaussian wavepacket 

vi>(r, 0) = i^—j exp{-a(r - qj^ + ,pT(r _ qj} (31) 

with qa,x = ^,(la,y = 0, Pa,x = 0,Pa,y = "1 and a = 0.25 is given in Fig. [1] Please 
note that we are not using an (anti-) symmetrized state here. Then both singlet and 
triplet states can be extracted from a single propagation without taking electron spin 
into account explicitly, see also [39] . 

We stress that without the Coulomb term the IVR results would be numerically 
exact and on top of the quantum ones (not shown). In the IVR case, with k = 0, 
even a single trajectory calculation according to the Thawed Gaussian approximation 
[33] is sufficient to generate the time-dependent semiclassical results. Therefore, the 
interesting case is the one with the Coulomb term. 

For K = 1, the results for some eigenvalues extracted from the time-series shown 
above are listed in Table [2] and are compared to the corresponding WKB results. We 
note that the IVR and WKB semiclassical results are coinciding to within the given 
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t 

Figure 1. Part of the auto-correlation time-series used to extract spectral information 
for 2d harmonium in a magnetic field (wl = ujq = 1) and with Coulomb interaction. 
Solid line: full quantum, dotted line IVR semiclassical result. 



m 


WKB 


IVR 


QM 





3.14 


3.14 


3.03 


1 


2.84 


2.84 


2.83 


2 


3.02 


3.02 


3.03 



Table 2. Comparison of WKB energies for = and different values of m with 
those from semiclassical IVR and full quantum calculations, including the Coulomb 
interaction and for wl = wq ~ 1. 

accuracy, although in the time-dependent calculations we did not(!) employ any Langer 
correction of the potential. For single electron dynamics in the Coulomb potential the 
Langer correction, necessary in the energy domain [30], is needed in the time-domain 
if one uses non- Cartesian coordinates [H] and seems to be unnecessary for Cartesian 
coordinates [121133], which is corroborated here. Furthermore, with better than 1 percent 
accuracy both semiclassical results are lying on top of the full quantum result, except for 
the m = eigenvalue. In the specific case considered, the full quantum result is around 3 
percent different from the semiclassical ones for m = 0. Furthermore, by coincidence, for 
the present parameters, two quantum eigenvalues for m = and m = 2 are degenerate. 
They do split up for a different choice of the ratio ul/ooq, however, where a similar 
agreement between the differently calculated energy eigenvalues is found (not shown). 
The marked difference (a few percent) between the full quantum result for m = and 
the semiclassical ones does persist. 
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5. Conclusions and Outlook 

We have shown that the time- dependent semiclassical IVR methodology of Herman and 
Kluk, based on Cartesian coordinates, is able to capture spectral features of two-electron 
quantum dots to a similar degree of accuracy as the WKB approach without the need to 
employ a Langer correction, complementing work on the comparison of WKB tunneling 
probabilities with Fourier transformed time-dependent semiclassical results (HJHS]. The 
semiclassical IVR method does not suffer from the restriction of WKB to integrable (e. 
g. one DOF) systems, however. This allows the investigation of non-circular quantum 
dots and also quantum dots in 3d or with more than two electrons. For the study of the 
dynamics of nuclear degrees of freedom the HK method, in the meantime, has become a 
widely used tool in the physics and chemistry communities (for an early review, see, e. 
g. [IS]). An open question in the application to more than two-electron systems is the 
problem of anti-symmetrization of the wavefunction, however. Furthermore, we stress 
that there is no restriction on the value of Q, as is the case for the exact analytical 
solutions of the problem presented by Taut |5], where the sequence of admissible Q 
starts with a value on the order of one and converges to zero. 
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